function lab8()
    % 读取图像并转换为灰度图像
    img1 = imread('lena.jpg');
    img2 = imread('photo.jpg');
    
    % 转换为灰度图像
    img1_gray = rgb2gray(img1);
    img2_gray = rgb2gray(img2);
    
    % 内置canny边缘检测
    img1_canny = edge(img1_gray, 'canny');
    img1_mycanny = mycanny(img1_gray, 0.03, 0.7, 0.3);
    
    img2_canny = edge(img2_gray, 'canny');
    img2_mycanny = mycanny(img2_gray, 0.07, 0.9, 0.2);

    % 显示图像
    set(0, 'defaultFigurePosition', [0, 200, 1500, 500]);
    figure;
    subplot(1, 3, 1); imshow(img1); title("原始图像");
    subplot(1, 3, 2); imshow(img1_canny); title("内置canny边缘检测");
    subplot(1, 3, 3); imshow(img1_mycanny); title("自行编写的canny边缘检测");
    
    set(0, 'defaultFigurePosition', [50, 150, 1500, 500]);
    figure;
    subplot(1, 3, 1); imshow(img2); title("原始图像");
    subplot(1, 3, 2); imshow(img2_canny); title("内置canny边缘检测");
    subplot(1, 3, 3); imshow(img2_mycanny); title("自行编写的canny边缘检测");

    % 保存结果
    imwrite(img1_mycanny, "edge1.tif");
    imwrite(img2_mycanny, "edge2.tif");
end

% 输入的参数分别是图像，低阈值，高阈值，二值化阈值
function f = mycanny(img, l, h, k)
    img_double = im2double(img);
    [width, height] = size(img);

    % 高斯模糊
    gauss = [1, 2, 1; 2, 4, 2; 1, 2, 1] / 16;
    img_blur = conv2(img_double, gauss, 'same');  

    % Sobel算子求梯度
    sobelx = [-1, 0, 1; -2, 0, 2; -1, 0, 1];
    sobely = sobelx';
    gradx = conv2(img_blur, sobelx, 'same');
    grady = conv2(img_blur, sobely, 'same');  
    M = abs(gradx) + abs(grady);
    theta = atan(grady ./ gradx);

    % 非极大值抑制
    N = zeros(size(M));
    for i = 2:width-1
        for j = 2:height-1
            dirc = theta(i,j);

            if abs(dirc) <= pi/8
                if M(i,j) == max([M(i,j-1), M(i,j), M(i,j+1)])
                    N(i,j) = M(i,j);
                end
            elseif abs(dirc) >= 3/8*pi
                if M(i,j) == max([M(i-1,j), M(i,j), M(i+1,j)])
                    N(i,j) = M(i,j);
                end
            elseif dirc > pi/8 && dirc < 3/8*pi
                if M(i,j) == max([M(i-1,j-1), M(i,j), M(i+1,j+1)])
                    N(i,j) = M(i,j);
                end
            elseif dirc > -3/8*pi && dirc < -pi/8
                if M(i,j) == max([M(i+1,j-1), M(i,j), M(i-1,j+1)])
                    N(i,j) = M(i,j);
                end
            end
        end
    end

    % 双阈值检测和边缘连接
    f = zeros(width, height);
    TH = h * max(max(N));
    TL = l * max(max(N));
    for i = 2:width
        for j = 2:height
            if N(i,j) > TL && N(i,j) < TH
                if N(i,j) > k   % 二值化
                    f(i,j) = 1;
                end
            end
        end
    end
end
